---
title: "Replication of analyses reported in Montero-Melis & Jaeger (2019) in *BLC*"
author: '[Guillermo Montero-Melis](https://www.mpi.nl/people/montero-melis-guillermo)'
date: '`r as.character(format(Sys.Date(), format="%d/%m/%Y"))`'
output:
  html_document:
    depth: 2
    number_sections: yes
    theme: default
    toc: yes
---


Introduction
============

This knitr document accompanies the paper

Montero-Melis, G., & Jaeger, T. F. (2019). Changing expectations mediate
adaptation in L2 production. Submitted to *Bilingualism: Language and Cognition*

For details about the paper, please refer to the main text and the Supplementary
Online  Material.
For details about how this report was generated, see the README file.


Set up working space
====================

Load libraries
--------------

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, fig.height = 4, fig.width = 5)
```

```{r, include = TRUE, echo = TRUE, warning = FALSE, message = FALSE}
library(dplyr)
packageVersion('dplyr')
library(tidyr)  # useful for gather() to convert from wide to long
packageVersion('tidyr')
library(mgcv)  # GAMs and GAMMs (Wood 2006)
packageVersion('mgcv')
library(itsadug)
packageVersion('itsadug')  # Visualizing GAMMs
library(lme4)
packageVersion('lme4')
library(ggplot2)
packageVersion('ggplot2')
library(GGally)  # for ggpairs()
packageVersion('GGally')
library(arm)  # used in script "plot_glmm_fnc.R"
packageVersion('arm')
library(knitr)  # for kable()
packageVersion('knitr')
```


```{r}
# Make sure you are in the folder where all the scripts and data are located,
# otherwise the script won't run.
getwd()
```


Convenience functions
---------------------

Because some models take a long time to fit, we use the self-mande convenience
function `load_or_fit`. This function first checks whether a model with a given
name already exists as a `.rda` file in the subfolder `fitted_models/`. If it does,
it loads it, otherwise it fits it and saves it to that folder.

```{r}
# Source function used to load models if they have already been saved, rather
# than fitting them anew each time the script is called.
source("load_or_fit_fnc.R")
# This function expects a subfolder "fitted_models/"; create it
dir.create(file.path("fitted_models/"))
# Create subfolder to save figures
dir.create(file.path("myfigures/"))
```

Format table for GLMM output

```{r}
# create a neat table of the summary of fixed effects of a mixed model	
glmm_tb <- function(fm) {	
  m <- round(summary(fm)$coefficients, 3)	
  tb <- as.data.frame(m)	
  names(tb) <- c("Estimate", "SE", "z-value", "p-value")	
  kable(tb)	
}	
```


Plot GLMM estimates (baseline condition):

```{r}
# Function that plots mixed model-estimates in logit space from baseline
# conditions, including speaker estimates
source("plot_glmm_fnc.R")
```


Plot GAMM estimates:

```{r}
# Two functions to a) plot adaptation effects for native and L2 speakers from
# GAMMs, and b) plot the effects by L2 speakers' proficiency.
source("plot_gams_fnc.R")
```


Global variables for figures
----------------------------

```{r}
# adjust figure heght/width when not going with default (espec. for 2x2 plots)
myfighe_NS_L2 <- 6
myfighe_L2_prof <- 6
myfigwi <- 7
```


Load data and set up data frames for analyses
-------------------------------------------

```{r}
# load data
d <- read.csv('data.csv', fileEncoding = 'UTF-8', stringsAsFactors = TRUE)
# Make sure Subject is a factor
d$Subject <- factor(d$Subject)
# in Group, let the native speakers (NS) be the reference group
d$Group <- factor(d$Group, levels = c('NS', 'L2'))
head(d, 3)
# load participant information 
# NB: Accidentally, no audio data was recorded for Subject 14 (L2), so this
# participant is omitted from all participant descriptors as well.
ppts <- read.csv("participant-info.csv", fileEncoding = "UTF-8", 
                 stringsAsFactors = TRUE)
head(ppts, 3)
tail(ppts, 3)  # L2 proficiency-related variables only available for L2ers
```

### Set up data for models (GLMMs and GAMMs)

We will do two things:

1) In each model, we want to allow for two comparisons: a) Path verb use in the
Baseline condition vs Path verb use in the Path-primed condition, and b) Manner
verb use in the Baseline condition vs Manner verb use in the Manner-primed
condition.
2) In addition, we want to define the correct data type and coding schemes for the
different factors.

```{r}
# Convert data to long format: each description will now have two rows, one for
# each value of VerbType (P_V = Path verb; M_V = Manner verb); note that these
# two rows were previously two different columns showing whether what was 
# encoded in the main verb in a given description encoded was path, manner,
# both or neither (see data frame d above).
d_long <- gather(d, VerbType, Used, P_V : M_V)

# Create copy of long data set, which we will further process for model fitting
d_mod <- d_long

# Combine VerbType (Path/Manner) and Condition (baseline/primed) into a single
# factor -- this is needed for the GAMMs. (The way to analyze crossed factors
# is to combine their levels into a single factor, see Wood's comment here:
# http://grokbase.com/t/r/r-help/113qaadxt4/r-how-to-add-in-interaction-terms-in-gamm)
d_mod$VbType_Cond <- with(d_mod, interaction(VerbType, Condition))	

# In order to make the relevant comparisons (see point 1 here above), baseline
# participants have to be compared twice, once for their path verb use to path-
# primed participants, once for their manner verb use to manner-primed 
# participants. This means that each of the two rows defined by VerbType will
# be used in one of the comparisons. The following renaming allows us to
# fit random by-subject smooths for each of those comparisons:
d_mod$Subject <- with(d_mod, interaction(Subject, VerbType))	

# Subset data for model fitting:
# We remove observations that	contain irrelevant comparisons given our research
# question, namely those that correspond to path verbs produced in the manner-
# primed condition or to manner	verbs produced in the path-primed condition.
d_mod <- d_mod %>% dplyr::filter(! VbType_Cond %in% c("P_V.Manner", "M_V.Path"))
# drop unused factor levels and order levels correctly
d_mod$VbType_Cond <- factor(d_mod$VbType_Cond, 
                            levels = c("P_V.Baseline", "P_V.Path", "M_V.Baseline",
                                       "M_V.Manner"))
# drop unused factors for subject (important for random effects estimation)
d_mod$Subject <- factor(d_mod$Subject)
# Condition recoded as a binary factor (Path/Manner become "Primed")
d_mod$Condition_bin <- d_mod$Condition
levels(d_mod$Condition_bin)[levels(d_mod$Condition_bin) %in% c("Path", "Manner")] <- "Primed"	
```


Coding scheme for factors:

```{r}
# define factor coding scheme (contrast coding)
# group
contrasts(d_mod$Group) <- - contr.sum(2)
colnames(contrasts(d_mod$Group)) <- "L2_vs_NS"
contrasts(d_mod$Group)
# Verb type
d_mod$VerbType <- factor(d_mod$VerbType)
contrasts(d_mod$VerbType) <- contr.sum(2)
colnames(contrasts(d_mod$VerbType)) <- "M_vs_P"
contrasts(d_mod$VerbType)
# Condition_bin contrast coding	
contrasts(d_mod$Condition_bin) <- - contr.sum(2)
colnames(contrasts(d_mod$Condition_bin)) <- "Primed_vs_Baseline"	
contrasts(d_mod$Condition_bin)	
```

Show data frame

```{r}
head(d_mod)
str(d_mod)
```



Participant descriptors
=======================

## As reported in the main text (Method section)

Number of participants by group:

```{r}
# Number of participants by group that go into the analysis (note data for
# one L2 speaker is missing because of recording failure)
table(ppts$Group)
```


Participant age:

```{r}
# age
ppts %>%
  group_by(Group) %>%
  summarise(Mage = mean(Age, na.rm = T), SDage = sd(Age, na.rm = T))
```


Learners' age of onset for learning Spanish 

```{r}
# Mean and SD
ppts %>%
  dplyr::filter(Group == "L2") %>%
  summarise(M  = round(mean(L2_AoO, na.rm = TRUE), 2), 
            SD = round(sd(L2_AoO, na.rm = TRUE), 2))
# range
range(ppts$L2_AoO, na.rm = TRUE)
```



Scores on offline cloze test ('proficiency score'):

```{r}
# Cloze scores
round(mean(ppts$ClozeScore, na.rm=T), 1)
round(sd(ppts$ClozeScore, na.rm=T), 1)
# By condition
ppts %>% 
  dplyr::filter(Group == "L2") %>% 
  group_by(Condition) %>%
  summarise(M_profic = mean(ClozeScore), 
            SD = sd(ClozeScore))
```


Difference in L2 proficiency across conditions (one-way ANOVA):

```{r}
# one-way ANOVA
ppts %>% 
  dplyr::filter(Group == "L2") %>%
  aov(ClozeScore ~ Condition, data = .) %>%
  summary
```


Computerized vocabulary task, which was deliberately easy (see main text).
(Data for this task is missing from one participant due to experimenter error.)

```{r}
# Scores on vocabulary task
round(mean(ppts$VocabTaskAcc, na.rm=T), 3)
round(sd(ppts$VocabTaskAcc, na.rm=T), 4)
# By condition
ppts %>% 
  dplyr::filter(Group == "L2") %>% 
  group_by(Condition) %>%
  summarise(M_profic = mean(VocabTaskAcc, na.rm = TRUE), 
            SD = sd(VocabTaskAcc, na.rm = TRUE))
```



## Participant information reported in the Supplementary Online Information (S2)

### Measures of L2 proficiency


Distribution of Cloze scores

```{r}
ppts %>% dplyr::filter(Group == "L2") %>%
  ggplot(aes(x = ClozeScore, colour = Condition)) +
  geom_density() +
  geom_rug() +
  facet_grid(Condition ~ .) +
  theme_bw()
```


Numerical summaries of L2 variables potentially related to L2 proficiency:

```{r, warning=FALSE}
ppts %>% dplyr::filter(Group == "L2") %>% 
  dplyr::select(ClozeScore, L2_SelfRating, L2_AoO, L2_instruction_years, VocabTaskAcc) %>%
  summary
```


Correlation between different L2-related measures that were collected:

```{r, warning=FALSE, fig.width=9, fig.height=9}
my_dens <- function(data, mapping, ...) {
  ggplot(data = data, mapping=mapping) +
    geom_density(..., alpha = 0.7) 
}
# Improved correlation matrix with ggpairs
ppts %>% 
  dplyr::filter(Group == "L2") %>% 
  dplyr::select(ClozeScore, L2_SelfRating, L2_AoO, L2_instruction_years, VocabTaskAcc,
                Condition) %>%
  GGally::ggpairs(
    mapping=ggplot2::aes(color = Condition),
    diag = list(continuous = my_dens),
    upper = list(continuous = wrap('cor', method = "spearman", size = 2.75, hjust=0.8)),
    lower=list(continuous="smooth")
    ) +
  # removes grid from correlations but not from other panels:
  theme(panel.grid.major = element_blank())
```


Do L2 participants in the different conditions differ along any of these measures?


```{r}
## Use Kruskal-Wallis test for this data 
# Cloze score (again)
ppts %>% 
  dplyr::filter(Group == "L2") %>%
  kruskal.test(ClozeScore ~ Condition, data = .)
# L2 self-rated proficiency
ppts %>% 
  dplyr::filter(Group == "L2") %>%
  kruskal.test(L2_SelfRating ~ Condition, data = .)
# L2 Age of onset
ppts %>% 
  dplyr::filter(Group == "L2") %>%
  kruskal.test(L2_AoO ~ Condition, data = .)
# Self-reported years of formal instruction
ppts %>% 
  dplyr::filter(Group == "L2") %>%
  kruskal.test(L2_instruction_years ~ Condition, data = .)
# Accuracy in vocabulary task
ppts %>% 
  dplyr::filter(Group == "L2") %>%
  kruskal.test(VocabTaskAcc ~ Condition, data = .)
```



### Other background information

Gender distribution of the participants:

```{r}
with(ppts, table(Group, Gender))
```


Age distribution:

```{r, fig.width=6, warning=FALSE}
ppts %>% 
  ggplot(aes(x = Age, colour = Condition)) +
  geom_density() +
  geom_rug() +
  facet_grid(Condition ~ Group) + 
  theme_bw()
```


Results and discussion
======================


Baseline condition: Native and L2 speakers’ lexical preferences in the absence of priming
------------------------------------------------------------------------------------


### Set up data set for GLMM


```{r, echo = T}
# Subset data from baseline condition only
d_basel <- d_long %>% 
  dplyr::filter(Condition == "Baseline") %>%
  dplyr::select(Subject : ClozeScore, VerbType, Used)
# Use contrast coding to compare groups...
contrasts(d_basel$Group) <- - contr.sum(2)
colnames(contrasts(d_basel$Group)) <- "L2_vs_NS"
contrasts(d_basel$Group)
# ... and Verb type
d_basel$VerbType <- factor(d_basel$VerbType)
contrasts(d_basel$VerbType) <- contr.sum(2)
colnames(contrasts(d_basel$VerbType)) <- "M_vs_P"
contrasts(d_basel$VerbType)
# show
head(d_basel, 3); tail(d_basel, 3)
```



### Reported in main text

Table 3 in main text:

```{r}
# Descriptive data for all analysed conditions (only relevant comparisons).
# Note the table was appropriately formatted in report.
d_mod %>%
  group_by(VerbType, Group, Condition) %>%
  summarise(Occurrences = sum(Used),
            TotalN = n(),
            Percentage = round(100 * sum(Used) / n(), 1)) %>%
  kable
```



```{r, echo = TRUE}
# load model or fit
# First define the expression to be passed to load_or_fit()
glmm_basel.expr <- 
"glmer(Used ~ VerbType * Group + (1 + VerbType | Subject) +
  (1 + VerbType * Group | VideoName),
data = d_basel, family = 'binomial',
control = glmerControl(optimizer='bobyqa', optCtrl=list(maxfun=2e5)))"
# Now pass it to the function (see source code for the inner workings)
load_or_fit("glmm_basel", glmm_basel.expr)
```


Fixed effects estimates:

```{r}
# output table
glmm_tb(glmm_basel)
```


Plot the results:

```{r}
# the function is defined in analysis/functions/plot_glmm_fnc.R
plot_basel(glmm_basel, d_basel, nb_sims = 500)
```

```{r}
# save to disk
ggsave("myfigures/Figure3.pdf", height = 4, width = 5)
ggsave("myfigures/Figure3.tiff", height = 4, width = 5, units = 'in', dpi = 600)
```





### Reported in Supplementary Online Material

Model summary:

```{r}
summary(glmm_basel)
```



Trial-by-trial adaptation of native speakers (Question 1)
---------------------------------------------------------

### Set up data frame:

```{r}
# Subset data to native speaker data only
d_ns <- d_mod %>% dplyr::filter(Group == "NS")
# drop unused factors for subject (important for gam fitting)
d_ns$Subject <- factor(d_ns$Subject)
# note log-odds of path verbs in baseline condition becomes the reference level
contrasts(d_ns$VbType_Cond)
```


### GAMM specification and fitting


First, we fit a more complex GAMM with by-subject and by-item factor smooths
for Trial (`s(Trial, Subject, bs = 'fs')` and `s(Trial, VideoName, bs = 'fs')`,
respectively).

```{r}
# the expression that is passed to load_or_fit()
gam_ns_itemsmooth.expr <- 
"bam(Used ~ VbType_Cond + s(Trial, by = VbType_Cond) + 
  s(Trial, Subject, bs = 'fs') + s(Trial, VideoName, bs = 'fs'),
data = d_ns,
family = 'binomial')"
# load model or fit (fitting took about 10 minutes on my laptop)
load_or_fit("gam_ns_itemsmooth", gam_ns_itemsmooth.expr)
```

Model summary:

```{r}
summary(gam_ns_itemsmooth)
```

We can see that the by-subject factor smooth for Trial (`s(Trial, Subject)`)
captures a fair amount of variance, whereas the by-item factor smooth for Trial
(`s(Trial, VideoName)`) does not.
Considering the interpretation of each of these factor smooths, this makes sense.
By-subject factor smooths capture how individual speakers are modifying
their encoding choices as the experiment progresses. One would indeed expect
those preferences to change at the individual level above and beyond the
population-level estimates. For the items, however, the factor smooth is modelling
the same type of changes at the level of the items (i.e., the individual events
being seen). Intuitively one expects less variation here: one would expect that
some items elicit more path encoding while others more manner encoding, but it
is less likely that the encoding preference associated to an item *changes* over
the course of the experiment. That effect should a priori be small if at all
existent. This is how we interpret the small variance associated to the by-item
factor smooth for Trial.

More practically, since factor smooths will take up many degrees of freedom and
slow down model fitting (here, it took twice the amount of time), let us compare
this model with another that has by-item random intercepts instead:

```{r gam_ns}
# the expression that is passed to load_or_fit()
gam_ns.expr <- 
"bam(Used ~ VbType_Cond + s(Trial, by = VbType_Cond) + s(Trial, Subject, bs = 'fs') +
  s(VideoName, bs = 're'),
data = d_ns,
family = 'binomial')"
# load model or fit (fitting took about 4.5 minutes on my laptop)
load_or_fit("gam_ns", gam_ns.expr)
```


```{r}
summary(gam_ns)
```

Note that the two models are identical in all coefficient estimates, except for
a very subtle difference for the estimated random effects of items (last line:
`VideoName`). We thus keep the simpler model.


In summary, the model is specified as:

`DV ~ VbType_Cond + s(Trial, by = VbType_Cond) + s(Trial, Subject, bs = 'fs') + s(VideoName, bs = 're')`

Predictors are:

- `VbType_Cond`: VerbType-Condition interaction as a fixed effect. The four levels
of this variable decompose into two crossed factors:
    - `VerbType`: Indicator variable that defines what is being measured by the DV,
  either the occurrence of path verbs or of manner verbs
    - `Condition`: Primed vs baseline
- `s(Trial, by = VbType_Cond)`: A smooth function of `Trial` allowing the function to differ for each level of Group-Condition (thin plate regression splines, the default)
- `s(Trial, Subject, bs = 'fs')`: Factor smooths for `Subject` to capture non-linear random effects of speakers
- `s(VideoName, bs = 're')`: Random intercepts by items (i.e., `VideoName`)



### Reported in main text


Plot differences between conditions (model estimates):

```{r, fig.height = myfighe_NS_L2, fig.width = myfigwi, echo = TRUE, results = 'hide'}
plot_gam_main(gam_ns, "NS")
# On my old computer, the above call threw an warning/error; this was solved by
# the following argument setting:
# plot_gam_main(gam_ns, "NS", mark.diff = FALSE)
```

```{r}
# Save to disk
tiff("myfigures/Figure4.tiff", width = 6, height = 5, units = "in", res = 600)
plot_gam_main(gam_ns, "NS")
dev.off()
pdf("myfigures/Figure4.pdf", width = 6, height = 5)
plot_gam_main(gam_ns, "NS")
dev.off()
```



### Reported in Supplementary Online Material

#### GAMM output

Summary of the GAMM for trial-by-trial adaptation in native speakers:

```{r, results = "asis"}
itsadug::gamtabs(gam_ns, type = "HTML")
```


#### Follow-up analysis (GLMM)

Rationale:

Based on the GAMM analysis, we can see that for native speakers there was a
significant adaptation effect for Manner, but not for Path verbs. However, 
we cannot conclude from one pattern being significant and the other not that
the two *differ* significantly. 
Unfortunately, this type of interaction analysis in a factorial design is 
difficult to implement in the GAM framework we are using (see, e.g.,
[here](http://grokbase.com/t/r/r-help/113qaadxt4/r-how-to-add-in-interaction-terms-in-gamm)).
We therefore carry out a follow-up mixed logistic regression analysis to address
this interaction, while noting that this is not an ideal solution:
we lose power, either because we model the effect wrongly (assuming linerarity
in a relation that is non-linear) and thus we may not capture the actual signal,
or because we reduce the number of observations to include only data points
that *do* have a linear relationship between predictors and dependent variable.
Here, we adopt this latter approach.

Since GLMMs assume linearity in log-odds space, We have to subset the data
to use only those trials that show an approximately linear relation with the
outcome variable (based on visual inspection). We keep trials 15 through 32.

Subset data:

```{r}
d_ns_lin <- d_ns %>% dplyr::filter(Trial >= 15)
```


Verify factors are contrast coded and `cTrial` is centred:

```{r}
# Verify factors are contrast coded:
contrasts(d_ns_lin$Condition_bin)
contrasts(d_ns_lin$VerbType)
d_ns_lin$cTrial <- d_ns_lin$Trial - mean(d_ns_lin$Trial)  # centre
round(mean(d_ns_lin$cTrial), 7)  # check
```


In these and all the following GLMMs, we try to stay conceptually close to the
random effects structure used in the GAMMs. This means we will generally
try to fit a model with random by-item and by-subject intercepts, as well as
a random by-subject slope for `cTrial` (corresponding to the smooth term
`s(Trial, Subject, bs = 'fs')` in the GAMMs).

Fit the model:

```{r}
# Random by-subject/item intercepts and random by-subject slope for cTrial.
# Expression to be passed to the load_or_fit function:
glmm_ns_lin.expr <- 
"glmer(Used ~ Condition_bin * VerbType * cTrial + (1 + cTrial | Subject) +
  (1 | VideoName),
data = d_ns_lin, family = 'binomial',
control = glmerControl(optimizer='bobyqa', optCtrl=list(maxfun=2e5)))"
# load model or fit it
load_or_fit("glmm_ns_lin", glmm_ns_lin.expr)
```

Model summary:

```{r}
summary(glmm_ns_lin)
```





Trial-by-trial adaptation of L2 learners (Question 2)
-----------------------------------------------------

### Set up data frame:

```{r}
d_l2 <- d_mod %>% dplyr::filter(Group == "L2")
# drop unused factors for subject (I think this is important for gam fitting)
d_l2$Subject <- factor(d_l2$Subject)
# note log-odds of path verbs in baseline condition becomes the reference level
contrasts(d_l2$VbType_Cond)
```


### GAMM specification and fitting

We proceed analogously as for native speakers (Question 1). We start by fitting
a more complex random structure, we then inspect the coefficients, and if warranted,
we then simplify it and check if anything changes.

```{r}
# the expression that is passed to load_or_fit()
gam_l2_itemsmooth.expr <- 
"bam(Used ~ VbType_Cond + s(Trial, by = VbType_Cond) + 
  s(Trial, Subject, bs = 'fs') + s(Trial, VideoName, bs = 'fs'),
data = d_l2,
family = 'binomial')"
# load model or fit (fitting took about 13 minutes on my laptop)
load_or_fit("gam_l2_itemsmooth", gam_l2_itemsmooth.expr)
```

Model summary:

```{r}
summary(gam_l2_itemsmooth)
```


```{r gam_l2}
# the expression that is passed to load_or_fit()
gam_l2.expr <- 
"bam(Used ~ VbType_Cond + s(Trial, by = VbType_Cond) + s(Trial, Subject, bs = 'fs') +
  s(VideoName, bs = 're'),
data = d_l2,
family = 'binomial')"
# load model or fit (fitting took about 5 minutes on my laptop)
load_or_fit("gam_l2", gam_l2.expr)
```


```{r}
summary(gam_l2)
```

The comparison is very similar to what we observed with native speakers:
The two models are pretty much identical in all coefficient estimates, except
for a difference for the estimated random effects of items (last line: 
`VideoName`), which does not affect the rest of the model. We thus retain and
report the simpler model.


### Reported in main text

Plot differences between conditions (model estimates):

```{r, fig.height = myfighe_NS_L2, fig.width = myfigwi, echo = TRUE, results = 'hide'}
plot_gam_main(gam_l2, "L2")
```

```{r}
# Save to disk
tiff("myfigures/Figure5.tiff", width = 6, height = 5, units = "in", res = 600)
plot_gam_main(gam_l2, "L2")
dev.off()
pdf("myfigures/Figure5.pdf", width = 6, height = 5)
plot_gam_main(gam_l2, "L2")
dev.off()
```


### Reported in Supplementary Online Material

#### GAMM output

Summary of the GAMM for trial-by-trial adaptation in L2 speakers:

```{r, results = "asis"}
itsadug::gamtabs(gam_l2, type = "HTML")
```


#### Follow-up analysis

The rationale for this analysis is analogous to that for the follow-up
analysis on Question 1 above.

Since GLMMs assume linearity in log-odds space, We have to subset the data
to use only those trials that show an approximately linear relation with the
outcome variable (based on visual inspection). We keep trials 1 through 13.

Subset data:

```{r}
d_l2_lin <- d_l2 %>% dplyr::filter(Trial <= 13)
```

Verify factors are contrast coded and `cTrial` is centred:

```{r}
# Verify factors are contrast coded:
contrasts(d_l2_lin$Condition_bin)
contrasts(d_l2_lin$VerbType)
# recenter cTrial
d_l2_lin$cTrial <- d_l2_lin$Trial - mean(d_l2_lin$Trial) 
round(mean(d_l2_lin$cTrial), 7)  # check
```


Fit the model:

```{r}
# Random by-subject/item intercepts and random by-subject slope for cTrial.
# Expression to be passed to the load_or_fit function:
glmm_l2_lin.expr <- 
"glmer(Used ~ Condition_bin * VerbType * cTrial + (1 + cTrial | Subject) +
  (1 | VideoName),
data = d_l2_lin, family = 'binomial',
control = glmerControl(optimizer='bobyqa', optCtrl=list(maxfun=2e5)))"
# load model or fit it
load_or_fit("glmm_l2_lin", glmm_l2_lin.expr)
```


Model summary:

```{r}
summary(glmm_l2_lin)
```





Effect of L2 proficiency on trial-by-trial adaptation (Question 3)
------------------------------------------------------------------

### Data and processing

The data is the same as we used to address Question 2, only now we will use
the predictor `ClozeScore` to measure L2 proficiency.

First rows of the data set:

```{r}
head(d_l2)
```


### GAMM fitting

Model for L2 speakers including proficiency:

```{r}
# model with L2 proficiency
gam_l2prof.expr <- 
"bam(Used ~ VbType_Cond + te(Trial, ClozeScore, by = VbType_Cond) +
  s(Trial, Subject, bs = 'fs') + s(VideoName, bs = 're'),
data = d_l2,
family = 'binomial')"
# load model or fit (fitting took about 6.5 minutes on my laptop)
load_or_fit("gam_l2prof", gam_l2prof.expr)
```

Model summary:

```{r}
summary(gam_l2prof)
```


### Reported in main text

Plot differences between conditions (model estimates):

```{r, fig.height = 3, fig.width = 9, echo = TRUE, results = 'hide'}
plot_L2_profic_diff(gam_l2prof, primed_cond = "Path")
plot_L2_profic_diff(gam_l2prof, primed_cond = "Manner")
```

```{r}
# Save to disk
tiff("myfigures/Figure6a.tiff", width = 7, height = 2.3, units = "in", res = 600)
plot_L2_profic_diff(gam_l2prof, primed_cond = "Path")
dev.off()
tiff("myfigures/Figure6b.tiff", width = 7, height = 2.3, units = "in", res = 600)
plot_L2_profic_diff(gam_l2prof, primed_cond = "Manner")
dev.off()
pdf("myfigures/Figure6a.pdf", width = 7, height = 2.3)
plot_L2_profic_diff(gam_l2prof, primed_cond = "Path")
dev.off()
pdf("myfigures/Figure6b.pdf",width = 7, height = 2.3)
plot_L2_profic_diff(gam_l2prof, primed_cond = "Manner")
dev.off()
```



### Reported in Supplementary Online Material


Model output (properly formatted):

```{r, results = "asis"}
itsadug::gamtabs(gam_l2prof, type = "HTML")
```

Plot each comparison showing the two curves (model estimates):

```{r, fig.height = 4, fig.width = 9, echo = TRUE, results = 'hide'}
plot_L2_profic(gam_l2prof, primed_cond = 'Path')
plot_L2_profic(gam_l2prof, primed_cond = 'Manner')
```



Self-priming (Appendix S7)
=========================

One of the reviewers brought to our attention the issue of self-priming, 
in the following terms:

> Strictly speaking, each priming trial consists of two different types of
utterances: the written and orally repeated sentence, and the sentence
produced by the participant to describe the picture. I take it that in the
adaptation analysis, only the number of trials has been taken into
consideration (that is, each prime+target pair is counted as 1). Would the
result be any different if in addition to the number of primes, also the
number of path or manner constructions produced as descriptions by the
participant were taken into consideration as well?

To address the reviewer's question, we take the following approach:

1. We compute a new variable `SelfProduction` that counts, for each trial, the
(cumulative) number of times a speaker has produced the target verb type (path or
manner verb, depending on the value of the `VerbType` variable) in all trials
up to the last one. This variable is always zero on the first trial and
goes up to maximally 31 (in trial 32, the last trial) if the participant
has produced the target verb throughout the whole experiment.
2. We fit a GAMM identical to the ones reported in the main paper, except we
use `SelfProduction` as a predictor *instead* of `Trial`.

The logic is as follows:

If self-priming can account for the results we observed, then we should see
*no differences* between the groups in the analysis below. This is because,
when using `SelfProduction` as the predictor, speakers in the different
conditions are being treated the same: a speaker for whom, say, `SelfProduction`
= 7 in a given trial, has been self-primed seven times, irrespective of whether
he or she is in the baseline or in the primed condition.
Conversely, to the extent that we still see a difference between the conditions,
this difference cannot be due to self-priming, but has to be attributed to the
only factor that differed between conditions, namely the fact that primed
participants also had to read and repeat the priming sentences.



## Create new variable `SelfProduction`

Add participant's cumulative productions, `SelfProduction` as a predictor:

```{r}
# Add number of produced structures as predictor --------------------------

# The cumulative number of productions is equal to the number of times the target
# structure # was used up until the last trial. This is what the following function
# comptues, where v is the vector of a subject's uses of the structure ordered by trial:
my_cumsum <- function(v) {
  out <- c(0, cumsum(v)[ 1 : (length(v) - 1) ])
  out
}

# Add the following predictors:
# a) SelfProduction: number of times the structure has been used in just production 
# b) prod_compr: Trial + the number of times it has been used in production. 
d_mod <- d_mod %>%
  group_by(Subject) %>%
  mutate(SelfProduction = my_cumsum(Used))
head(data.frame(d_mod))
```


Q1 re-analysis
--------------

### Subset data

```{r}
# Subset data to native speaker data only
d_ns <- d_mod %>% dplyr::filter(Group == "NS")
# drop unused factors for subject (important for gam fitting)
d_ns$Subject <- factor(d_ns$Subject)
```


### Model with self-priming from own production only


```{r, results = 'hide'}
## Model with self-priming in production only

# the expression that is passed to load_or_fit()
gam_ns_prod.expr <- 
  "bam(Used ~ VbType_Cond + s(SelfProduction, by = VbType_Cond) + s(SelfProduction, Subject, bs = 'fs') +
s(VideoName, bs = 're'),
data = d_ns,
family = 'binomial')"
# load model or fit (fitting took about 4.5 minutes on my laptop)
load_or_fit("gam_ns_prod", gam_ns_prod.expr)
```

Model summary:

```{r}
summary(gam_ns_prod)
```

Plot:

```{r, fig.height = myfighe_NS_L2, fig.width = myfigwi, echo = TRUE, results = 'hide'}
plot_gam_main(gam_ns_prod, "NS", show_xaxis = 'SelfProduction')
```



Q2 re-analysis
--------------

### Subset data

```{r}
# Subset data to L2 speaker data only
d_l2 <- d_mod %>% dplyr::filter(Group == "L2")
# drop unused factors for subject (I think this is important for gam fitting)
d_l2$Subject <- factor(d_l2$Subject)
# note log-odds of path verbs in baseline condition becomes the reference level
contrasts(d_l2$VbType_Cond)
```


### Model with self-priming in production only

```{r, results = 'hide'}
## Model with self-priming in production only

# the expression that is passed to load_or_fit()
gam_l2_prod.expr <- 
  "bam(Used ~ VbType_Cond + s(SelfProduction, by = VbType_Cond) + s(SelfProduction, Subject, bs = 'fs') +
s(VideoName, bs = 're'),
data = d_l2,
family = 'binomial')"
# load model or fit (fitting took about 4.5 minutes on my laptop)
load_or_fit("gam_l2_prod", gam_l2_prod.expr)
```

Model summary:

```{r}
summary(gam_l2_prod)
```


Plot:

```{r, fig.height = myfighe_NS_L2, fig.width = myfigwi, echo = TRUE, results = 'hide'}
plot_gam_main(gam_l2_prod, "L2", show_xaxis = 'SelfProduction',
              ylim1 = c(-6, 6), ylim2 = c(-.5, 8))
```



Q3 re-analysis
--------------

We use the same data as in Q2.


### Model with self-priming in production only

```{r}
## Model with self-priming in production only

# the expression that is passed to load_or_fit()
gam_l2prof_prod.expr <- 
  "bam(Used ~ VbType_Cond + te(SelfProduction, ClozeScore, by = VbType_Cond) +
s(SelfProduction, Subject, bs = 'fs') + s(VideoName, bs = 're'),
data = d_l2,
family = 'binomial')"
# load model or fit (fitting took about X minutes on my laptop)
load_or_fit("gam_l2prof_prod", gam_l2prof_prod.expr)
```


Model summary:

```{r}
summary(gam_l2prof_prod)
```


Plot effects by condition:

```{r, fig.height = 4, fig.width = 9, echo = TRUE, results = 'hide'}
# Plot effects by condition
plot_L2_profic(gam_l2prof_prod, primed_cond = 'Path', show_xaxis = 'SelfProduction',
               ylim1 = c(-9, 8))
plot_L2_profic(gam_l2prof_prod, primed_cond = 'Manner', show_xaxis = 'SelfProduction',
               ylim1 = c(-9, 8))
```


Plot differences between conditions

```{r, fig.height = 3, fig.width = 9, echo = TRUE, results = 'hide'}
# Plot differences between conditions (model estimates):
plot_L2_profic_diff(gam_l2prof_prod, primed_cond = "Path", show_xaxis = 'SelfProduction',
                    ylim1 = c(-2.5, 11))
plot_L2_profic_diff(gam_l2prof_prod, primed_cond = "Manner", show_xaxis = 'SelfProduction',
                    ylim1 = c(-2.5, 11))
```





Session info
============

```{r}
sessionInfo()
```


